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Abstract In this paper, we consider an efficient iterative approach to the solution 
of the discrete Helmholtz equation with Dirichlet, Neumann and Sommerfeld-like 
boundary conditions based on a compact sixth order approximation scheme and pre- 
conditioned Krylov subspace methodology. A sixth order compact scheme for the 
3D Helmholtz equation with different boundary conditions is developed to reduce 
approximation and pollution errors, thereby softening the point-per-wavelength con- 
straint. The resulting systems of finite-difference equations are solved by different 
preconditioned Krylov subspace-based methods. In the majority of test problems, 
the preconditioned Generalized Minimal Residual (GMRES) method is the supe- 
rior choice, but in the case of sufficiently fine grids a simple stationary two-level 
algorithm proposed in this paper in combination with a lower order approximation 
preconditioner presents an efficient alternative to the GMRES method. In the anal- 
ysis of the lower order preconditioning developed here, we introduce the term "fc-th 
order preconditioned matrix" in addition to the commonly used "an optimal precon- 
ditioner". The necessity of the new criterion is justified by the fact that the condition 
number of the preconditioned matrix AA^ ' in some of our test problems improves 
with the decrease of the grid step size. In a simple ID case, we are able to prove 
this analytically. This new parameter could serve as a guide in the construction of 
new preconditioners. The lower order direct preconditioner used in our algorithms 
is based on a combination of the separation of variables technique and Fast Fourier 
Transform (FFT) type methods. The resulting numerical methods allow efficient im- 
plementation on parallel computers. Numerical results confirm the high efficiency of 
the proposed iterative approach. 
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1 Introduction 

In recent years, the problem of increasing the resolution of existing numerical solvers 
has become an urgent task in many areas of science and engineering. Most of the 
existing efficient solvers for structured matrices were developed for lower-order ap- 
proximations of partial differential equations. The need for improved accuracy of 
underlying algorithms leads to modified discretized systems and as a result to the 
modification of the numerical solvers (see e.g. 1241). 

The use of a lower order preconditioner for efficient implementation of high- 
order finite-difference and finite-element schemes has been under consideration for 
a long time (see e.g. lfT4l . IfTSl ). In this paper, a compact sixth order approximation 
finite-difference scheme is developed, and a lower order approximation direct solver 
as a preconditioner for an efficient implementation of this compact scheme for the 
Helmholtz equation in the Krylov subspace method framework is considered. This 
approach allows us to utilize the existing lower order approximation solvers which 
significantly simplifies the implementation process of the higher resolution numerical 
methods. 

The model problem considered in the paper is the numerical solution of the 
Helmholtz equation 

V^u + k^u = f, inQ, (I) 

with the Dirichlet, Neumann, and/or Sommerfeld-like (radiation) boundary condi- 
tions 

M = 0,on dQi , 

u„ = 0,on dQ2, (2) 
u„ — iku = 0, on 90,^, , 

where £2 = {Q< x^y^z < a}, A; is a complex valued constant, and (9i2i , d£22 and dQj, 
are different boundary sides of the rectangular computational domain Q. 

It is known that for a given error level in the numerical approximation to the solu- 
tion of the Helmholtz equation, the quantity {Re{k)Y^^ h'' needs to be constant, where 
p is the order of finite-difference scheme and h is the grid size. This phenomenon is 
known as "pollution" lUO. One way of reducing the pollution error is to increase the 
order of accuracy of the scheme. In this paper a sixth order compact finite-difference 
scheme is considered to address this problem. In the cases of the 3D Helmholtz and 
Dirichlet boundary conditions, this scheme was proposed by Sutmann in ETll . The 
2D version of this scheme with Dirichlet and Neumann boundary conditions was de- 
veloped in Ifr7ll22]| . In this paper, the method is extended to the explicit compact sixth 
order approximation of Neumann and Sommerfeld-like boundary conditions in the 
3D case. The extension of the known approach for the Neumann boundary condi- 
tions is straightforward. But in the case of the Sommerfeld-like boundary conditions 
the method is nontrivial and requires the introduction of a new auxiliary function. In 
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the case of the variable coefficient Helmholtz equation and Sommerfeld-type bound- 
ary conditions, the fourth order compact approximation scheme was considered in 
0. 

The resuhing discretization leads to a system of linear equations with block 27- 
diagonal structure. In general, the matrix of this system is neither positive definite 
nor Hermitian. Hence, most iterative methods either fail to converge or converge too 
slowly, which is impractical. For the solution of this system we propose using a com- 
bination of Krylov subspace-based methods and the FFT preconditioner Concerning 
other approaches for the solution of the problem ([T]!-©, we refer to Bayliss, Gold- 
stein and Turkel [2] for a preconditioned conjugate-gradient algorithm, Kim 1 16| for 
a domain decomposition method, and Douglas, Hensley and Roberts [6| for an ADI 
algorithm. A very efficient multigrid method based on a shifted-Laplace precondi- 
tioner for the Helmholtz equation is presented in |23|. The analysis of the multigrid 
algorithms in the case of nonsymmetric and indefinite problems can be found in |4|. 

On the other hand, the solution of this problem by a direct method based on 
Gaussian elimination requires a prohibitive amount of additional storage and com- 
puter time and thus has limited use. The most promising results in the solution of 
a similar problem have been obtained by preconditioned Krylov subspace methods 
in 191. In this paper we generalize some approaches developed for the second-order 
central difference discretization of the Helmholtz equation by Elman and O'Leary [3 
O and the author and others |12| to the case of compact sixth order approximation 
scheme. 

The key to the fast convergence of the suggested iterative method is the choice 
of the preconditioning matrix. In our algorithm we use a preconditioner based on 
the second order central difference approximation of the Helmholtz equation and 
corresponding boundary conditions. The inversion of the preconditioning matrix at 
each step of the Krylov subspace method is done by a direct solver based on the 
FFT technique which requires 0{N^ logA^) operations, where is the number of grid 
points in each direction. 

Numerical experiments with test problems demonstrate the high resolution of the 
sixth order compact finite-difference scheme, as well as the computational efficiency 
of our preconditioned Krylov subspace type numerical method. In most situations, 
the GMRES method demonstrates the best convergence properties. However, in the 
case of sufficiently fine grids, we propose using a simple stationary two-level method 
(see e.g.pO]). This method was naturally constructed in the analysis of the conver- 
gence of the GMRES algorithm. So, the choice of the parameters in this method is 
not based on the minimization of the spectral radius of the iteration matrix on each 
step but rather on the construction of a particular linear combination from a Krylov 
subspace. This is why, for the purpose of this paper, we took the liberty of calling this 
approach the Simplified Krylov Subspace (SKS) method. On sufficiently fine grids 
this method requires much less processor time than the GMRES method, though the 
number of iterations until convergence is still larger than in the case of the GMRES 
method as should be expected. Also, since the implementation of SKS algorithm does 
not require the calculation of scalar products, it has greater potential than the GMRES 
method for implementation on parallel computers. A distinguishing feature of these 
approaches is that the number of iterations required for the convergence of Krylov 
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subspace iterations decreases as the size of the discretized system increases. To ex- 
plain this fact we introduce "the order of the preconditioned matrix" as a parameter 
to quantify the rate of convergence of the condition number of the preconditioned 
matrix AAp^ to 1 on a sequence of grids. In some simple situation, we were able to 
find this parameter analytically. We believe that this parameter is more informative 
then the commonly used "an optimal preconditioned' (see e.g. ifTSll . p. 196) in the case 
of a lower order preconditioners and may be used as guide in the further development 
of preconditioners of similar type. But we must notice that even in this paper this 
parameter has limited application. 

This conclusion is based on the theoretical analysis of some simple situations 
and confirmed in our numerical experiments. Some preliminary numerical results on 
the application of this approach were presented at the 10th International Conference 
on Mathematical and Numerical Aspects of Waves, Vancouver, Canada, 24-29 July, 
2011 ifTTl . 

The rest of the paper is organized as follows. In Section 2, the main idea of the 
proposed sixth order approximation compact method is presented in the case of the 
ID problem. To analyze the convergence of the algorithms developed here, variations 
of known convergence estimates for the Krylov subspace methods are considered. In 
this section, simplified approaches based on the Krylov subspace iterations are also 
presented. Section 3 focuses on the development of the sixth order compact approx- 
imation scheme in the case of Neumann and Sommerfeld-like boundary conditions. 
In Section 4, the effectiveness of the proposed algorithms is demonstrated on a series 
of test problems. 



2 A one-dimensional model problem 

Let A and Ap be matrices derived from six and second order approximations to ([T]i and 
the Dirichlet boundary condition (|2| using a meshx,- — ih,i = 0,1,..., N+ I, li^ TTfT' 
The standard notation for the first and second order central differences at ith grid 
point is given by 

Mi+l— M(-l e2 M,+ l — 2m, + M,_1 

' ^ 2h ' ^ Jfi ' 

where m, = m(x,). The difference operators 5y, 5^, 5^ and 5^ used in the following 
sections are defined similarly. The sixth order approximation to the second derivative 
at the ith grid point can be written as 

«;' = 5^-^«r)-^«r)+0(/,^). (4) 

As usual, in the case of compact schemes, we need to use the original equation to find 
appropriate relations to eliminate the fourth and sixth derivatives in (|4|. By using the 
second order central difference of the fourth derivative of m at jc = x, 



(5) 
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and the expression for the fourth derivative of u{x) from ([T]i 

u^^^-k'u'f+f;', (6) 

the relation (|4|i can be expressed in the form 

1 - ^(1 + ^5^)) «r = 5^u, - ^ (l + ^5,^) + 0{h% (7) 

where// = f{xi). After using the discretized system corresponding to the compact 
sixth order approximation scheme for ([T])-© can be presented as 

diUi-x+diUi + d^Ui+i ^Fi, 1,...,A^, 
t/o = 0, Un+x = 0, 

VA--*/!-* \ f k'^h* If If \ I Ih'* fit I ( fll 



where f/, is the sixth order discrete approximation to m(x,). This system can also be 
rewritten in the form AU = F, where 



A = dxA +d2l 



(8) 



and 



A = 



1 

1 1 







1 1 
1 



(9) 



In a similar way the preconditioning matrix based on the second order central differ- 
ence approximation can be presented as 



Ap=A-(2-r/!^)/. 
Finally the right preconditioned system can be written 



AA-'Y=F, 
ApU = Y. 



(10) 



(11) 



The main goal of this paper is to demonstrate the computational efficiency of a 
preconditioning technique based on the use of a lower order approximation discrete 
system in the numerical implementation of higher order compact finite-difference 
schemes. To consider this construction in more general form, we introduce the defi- 
nition of a kth order preconditioned system. 
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Definition 1 Given two nonsingular x matrices A and Ap, a system in the form of 
(fTTT i is said to be a kth order preconditioned system if the NxN matrix AA^ ' is diago- 
nahzable and can be expressed as V^^{I + h'^D)V, where h <1,D =diag{dii, ...,d„„) 
with max|t/„| < M, where M does not depend on h; and V is an x matrix of 

eigenvectors of AAp ' . 

Using this definition we present the convergence analysis of the GMRES method 
appHed to the system (fTTT i in the following theorem. 

Tlieorem 1 Let a system in the form of ( 1771 ) be a kth order preconditioned system. 
Then the nth iteration t/^"' of the GMRES method applied to this system satisfies the 
convergence estimate 

|k„||2<K-2(V)(M/!^)"|ko||2, (12) 

where KiiV) = HV"' Ibll V||2 and r„ ^ F - AU^"l 

Proof. Since the matrix AA^ ' is diagonalizable, it is well known (see e.g. llT0l . lfT9l ') 
that the residual of nth iteration of GMRES algorithm satisfies 

lk„||2< K-2(y) min max \p{l +h%i)\-\\ro\\2. (13) 

peP„.,p(0) = l i=\....,N 

n-l 

Consider the polynomial of nf/; degree Pn{x) ~ 1 — X3'„_i(x), where i (x) = ^ Oc['_|_jX-' 

*:=0 

and the coefficients satisfy the equation Pn{\ +x) = —a"x",n > 1. It is easy to 
see that the coefficients a^', k—l,...,n satisfy 

L«r-i, (14) 

k=l 



and the unique solution of this system is given by = (—1)*^^' ('J^,k—l,...,n. Then 
the convergence estimation for the GMRES method becomes 

Ik«ll2 < K-2(V) min max |p(l+/;*^t/,,)|-|ko||2 (15) 
<K2{V) max \p„{^+h%)\-\\ro\\2 (16) 

i= 1 N 

< K2{V){Mh')"\\ro\\2. 

The proof is somehow trivial and could be significantly simplified but we present 
it in such a form to use later in the construction of the simplified iterative method 
based on this proof. 

In the same way we can derive another useful estimate expressed in the following 
corollary. 
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Corollary 1 Let the matrix AAp ' in the system ( 1771) be expressible as 

where D = diag{dii, ...,d„n), M = max|t/,-, | < 1 and V is the N x N matrix of eigen- 

i 

vectors ofAAp^. Then the nth iteration U„ of the preconditioned GMRES method 
applied to this system satisfies the convergence estimate 

lk„||2< K-2(V)M"||ro||2. (17) 

We can derive another well known estimate based on a standard application of Cheby- 
shev polynomials (see e.g. [19]). We consider only the case of real valued matrices. 
The following theorem provides the details of the proposed technique. 

Theorem 2 Let the matrix AA^^ in the system ( 1771) be diagonalizable and expressible 
in the form V^^ {I + D)V, where D = diag{dii, ...,d„„), —1 < m = min da < M = 

i 

max da and V is the N x N matrix of eigenvectors ofAA^K Then the nth iteration Un 

i 

of the preconditioned GMRES method applied to this system satisfies the convergence 
estimate 

Proof. This estimation immediately follows from the general complex valued ma- 
trix result (see e.g.09|, P-206). 

Next we will applied the developed convergence estimations to the analysis of the 
solution of the system (fTTT i. 



and min 

7=1 N 



Theorem 3 Let N xN matrices A and Ap be defined by dSf and dTOl l with h < 



It 



2 



fe2 



> 5o > 0, where 5q is constant. Then AAp is a second 



order preconditioned system and we have the convergence estimate. 



Proof. The matrices A and Ap have the same set {Vj,j = 1, ...,A^} of eigenvectors, as 
the matrix A from The eigenvectors are given by (see e.g.|8|) 

V] = V2hsm{jlKh)J,l = l,...,N. (20) 

The matrix V is unitary, so || V^' II2 = ll^lb ~ 1- The eigenvalues of the matrices Ap 
and A are given by 

= -4sin2 (^^^ +hV, for ./ = 1, ...,A? and (21) 

9 / iKh\ 99 A A (\A + cos (inh)) 

Xj = -4s n^ i — +hV-h^k^-^ — ^- — ^, for / = l,...,N 

^ \ 2 I 180 ' i ' ' 
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respectively. We can write the eigenvalues of the preconditioned matrix AA^ in the 
form 



3 >,4.4 (14+cos(,/;r/i)) 

^^1 + — m for j 

a; 4sin2(if^)-/z2fe2' 



1 A^. 



(22) 



The second term in the right hand side can be estimated using m < 



- 11 < M, 



where m = 13/^^/720 andM = h^k'^/iUdo). So AA„^ is a second order precondi- 
tioned matrix and using (fT2l) we obtain the estimate stated in the theorem. 

The first condition of this theorem hk < 2n/l0 is a common natural restriction 
that comes to play even when one wants to visualize a numerical solution for qual- 
itative analysis. It is just not an accurate representation of the solution if there are 
fewer then 5 points per half wave length. Moreover, to avoid the so-called the "pol- 
lution" phenomenon Bayliss et. al. [3J introduced the restriction that k''^^h'' should 
be constant, i.e. to maintain a fixed accuracy of a scheme, the number of grid points 
must grow as k^^P where p is the accuracy order of the scheme. So, if one satisfies 
the restriction that avoids "pollution" phenomenon, the first condition of the theorem 
is satisfied almost automatically and it is not much of a restriction at all. 

The second condition of the theorem is a requirement that avoids the discrete 
spectrum of the operator Since the resolvent set of a bounded linear operator is open, 
we can always do this by introducing a small change to k^ if necessary. 

Next, we present how we can choose 5o in some important cases. 

In the simplest case of k^ < 9, we can take 5o = 1/3. 

Consider how to get an estimate for 5q in the case of a given k > n and a se- 
quence of grids with the grid sizes hi > h2 > hi, . . . . The first restriction of the the- 
orem gives us the condition kh <jk. We can see that if the j in the argument of 



rmn 



4sin2(if^) 



/l2 



-k' 



were to take on values of all real numbers from 1 to A^, then 



the minimum of the expression under consideration would be zero. 

Now let's denote a{j) — ^ and consider for which a the expression is zero if 
kh = j^. We consider the worst-case scenario to justify this approach for all kh < '^^ 



It is trivial to find that Oq = sin 



-'(f) = sin-H^) 



10- 

.32. So, it is clear that the 



minimum of the function occurs either at jo = Li^J or ji — Jo + 1 since sin(a(j)) 
is increasing on the interval 1 < j < N. There exists a neighborhood of OCo which 
includes both a{jo) and a{ji) such that a < .64 which allows us to use the Taylor 
series to represent sin^(a) in this neighborhood. We use only first two terms in the 
series, i.e. sin^(a) — — |cos(2(^)a"^,0 < ^ < 0.64. Now we can substitute this 
expression in the equation 



mm 

j=i,...^ 



4sin2(M) 



mm 

j=joJi 



> min 



/;r2--cos(2^)(j;r)V-fe2 



k'^jt,n'Jin'-k'--ij7tr{h,y 



>5o. 
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So, we just need to avoid the values for k and h\ for which one of the terms in the 
brackets is nonpositive. In some situations we would need to choose a sufficiently 
small grid size step h\ for the coarsest grid in a sequence. For example, in the case 

;t = 20 and /!i = 1 /32, min ^400 - 6V,l^Tt^ - 400 - ^ ^) > min(44,63) = 5q. 

This ^3 is independent of h and can be used in our convergence analysis on a sequence 
of grids. 

It follows from ( fT9] l that the number of iterations for the Krylov subspace algo- 
rithm decreases as the grid size h decreases. This result proves that the algorithm 
developed here yields a very effective iterative technique for the implementation of 
a higher order approximation scheme. Later we will consider the extension of this 
method to 3D problems. 

The proof of Theorem [T] suggests a simplified version of the Krylov subspace 
method for the solution of (fTTl i. Indeed, using the expression for the exact solution 
of the system (fT4l i. we have a," = — 0!/,'^[ , for n= 1, ... and the following recurrence 
expression for y„ 

Jn (z) = 1 - zjn- 1 (z) + Jn- 1 (z) , whcre n > 1 . 

This allow us to derive the SKS algorithm described in the table labelled Algorithm 
1. 



Algorithm 1 Simplified Krylov Subspace (SKS) algorithm. 

1: Let F'"' = Ap't/'"' = be the initial approximation, let the tolerance be tol, and let the maximum 
number of iterations be M 



2 


ro = F,Y = ro and So = ||ro||2 


3 


y = i 


4 


while j <M and e < tol do 


5 




6 




7 


r = ro — w 


8 


5 = lkll2 


9 




10 


Y = Y + r 


11 


J = j+i 


12 


end while 


13 


C/'^' is the iterative solution of AC/ 



This is a simple stationary two-level method (see e.g.f24l). However, the choice 
of the parameters in this algorithm is not based on the minimization of the spectral 
radius of the iteration matrix on each step but rather on the construction of a particular 
linear combination of vectors from a Krylov subspace. This is why, for the purpose of 
this paper, we call this approach the Simplified Krylov Subspace (SKS) method. The 
proposed method does not require an orthogonalization procedure, i.e. requires no in- 
ner products, which is attractive for an implementation in a parallel computing envi- 
ronment. In addition, this method does not use estimates of maximum and minimum 
eigenvalues of the preconditioned matrix which is the drawback of the Chebyshev 
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acceleration algorithm (see e.g. ||T9| ). The Chebyshev acceleration method is natu- 
rally constructed in the proof of the Theorem |2] and we will compare the numerical 
effectiveness of the methods in Section 4. We should notice that in most numerical 
experiments, the GMRES method exhibits the best convergence properties, but in 
some situations the SKS algorithm proves to be more efficient. Some such examples 
are considered in Section 4. 

The SKS algorithm also gives a good criterion for evaluating the quality of a pre- 
conditioner in the form of the order of preconditioned matrix ( fTTI ). In this paper we 
will calculate this as follows. Consider two grids with the grid size h and 7/1, where 
7 > 1. Let the Z2^iiorm of the residuals of the SKS method on first two iterations 
be ||r[||2 and \\r2\\2 on the first grid and ||rj'''||2 and ||r2''||2 on the second grid. Let 
I -''II II ■''''I 

e'^ = jjY' ^iid ^ Ih ■ Now we can approximately calculate the order of pre- 
conditioned matrix by 

Ine'^'-lne'' 

We consider this parameter in the discussion of the results of the numerical experi- 
ments. 



3 Three dimensional problems 

3.1 A sixth order approximation compact scheme 

In this section we present a three dimensional compact sixth order approximation 
finite-difference scheme that was first introduced in [211, and we develop a sixth or- 
der compact explicit approximation of the Neumann and Sommerfeld-like boundary 
conditions (|2]i. In this discussion we consider a uniform grid i2/, = = 
ih,yj = jh,z = kh;i,j,k ^Q, ...,N - l,h ^ a/{N - I)}. 

First, we consider a finite-difference approximation of ^ in the form 



where Ujj^k = u{xi,yj,Zk) and = f{xi,yj,Zk) and 



(24) 



'^'J,k - 



d'^u d'^u d^u 
12 \ d7^^^J?^ 



360 \ d^^^^J^ 



0{h% 



ij,k 



Using the appropriate derivatives of ([TJ we can write the sixth order compact 
approximation of the Helmholtz equation in the form 



30 

k^h^ k^h* 



12 



Ui.jM = 



360 



(25) 
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1 + ]fiik + —[ 1 ^ fi i k + ^ fi i k - 

12 360 r 12 I 30/ ■"'^■'^ 360 



( d'f d^f d^f 
90 ydx^dy"^ dx^dz^ dy^dz^ 



For the two dimensional case this scheme was proposed in [ 17 |. In the three dimen- 
sional case with Dirichlet boundary conditions, it was developed in ll2n . We will 
consider the iterative implementation of this scheme based on preconditioned Krylov 
subspace type algorithms. 



3.2 Boundary conditions 

The implementation of the Dirichlet boundary conditions (O is straightforward but 
the explicit compact approximation of Neumann and Sommerfield-like boundary 
conditions requires careful consideration. The sixth order compact approximation 
of the Neumann boundary condition in the two dimensional case was considered in 
IfTTl . Here we extend this approach to the three dimensional problem. 



3.2.1 Neumann boundary conditions 

In this subsection the Neumann boundary conditions are considered in the form 



du 

Tz 



p{x,y). 



(26) 



We restrict our consideration to only one side of the computational domain Q. 
By using Taylor series, we can derive the sixth order approximation formula 



du 

^^"•Jfi = Yz 



i.J.O 



~6 a? 



d^U 

120 a? 



(27) 



ijfi 



To express the third and fifth derivatives in ( |27] | we can differentiate the original 
Helmholtz equation ([T]l, assuming sufficient smoothness of the solution and the right 
hand side. After substitution of these derivatives into dZTb . the resulting expression is 



+ 





du 




'df 




Tz 




_dz 




d^u 


d\ 



d\ 



d^u 



2 du 



dx^dz dy^dz dz 



120 

0{h^ 



dz^ 



dz^ dy? 



dz^dy^ 



d\ 



Ufi 

d\ 



dz dx^dz dy^dz 



^du 
" dz 



(28) 

ij,0 



Next, we approximate the third order mixed derivatives in the second bracket of ( |28] ) 
by a fourth order approximation formula and in the second line using a second order 
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approximation scheme, which yields 



1 - 

If- r 
6 



5^5^ II + 5t5„ m 





5m 


) 













;.;,o 



4- 



20 



dz 



( 



120 [ 



!.;,0 



^ T20 [^'^-^ "'•^■"^ ^ T20 

Now we can simplify this and write 



dx^dz^ dy^dz^ 



6 



20 



360 

7^ 
360 



6 120 



du 

Tz 



iJ.O 



20 



ijfi 



dzd^x dzd^y 







-,.^180 


dx^dy^dz 



dx^dz dy^dz 



iJ.Q 



Tso 



iJfi 



dx^dz dy'^dz 



120 dz^ 



0{h^). 



Finally, to complete the explicit scheme for the boundary conditions, we present the 
previous equation in the form 



'In 



30 



2;,2 



1 - 

360 
180 



h'-k 



6 



/2 V \ 
120 J 



5m 
5l 



;.;,o 



20 



^z 



120 



dz^ 



i.j.O 



dzd^x dzd^y 









dx^dy^dz 



dx^dz dy^dz 



ij.o 



180 



ij.o 



dx^dz dy^dz 



0{h'). 



Now by using (O and replacing m, y ^ with jt, we can write 



Ui.j.-i + ^il 



2/1 1 - ■ 



30 

4^4 



(29) 



+ 



Tso 



120 



30 

lflf\ \df 



dz 



60 



dz^ 



ijfl 



dzd^x dzd^y 



iJfl 



90 



+ ^ V — + — 

^jc^^y^ \ dx^ dy^ J dx^ dy^ 



iJfi 
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This formula allows us to eliminate the term Uij,-i in (IZST l. The explicit implemen- 
tation of the Neumann boundary conditions on the other boundary pieces can be 
conducted in a similar way. 



3.2.2 Sommerfield-like boundary conditions 

In this section we consider the implementation of the Sommerfeld-like boundary con- 
ditions (IS which are often used in scattering problems 



du 



' iku 



= 0. 



(30) 



:=0 



The difficulty in using similar methods to those described in the previous section can 
be seen in Equation ( |29] |. The last term in this equation requires calculation of the 
fourth derivative of the right hand side of ( l26b and if it depends on an unknown vari- 
able u, the use of a compact sixth order approximation scheme becomes problematic. 
To avoid such a difficulty, we introduce a new variable v ~ e'^^u. Then the Helmholtz 
equation becomes 



dv - 



where / = e'^'/. Now (|30] | can be rewritten in the form 



dz 



z=0 



The sixth order approximation of ( l32b becomes 



(31) 



(32) 



- Yz 



IJ.O 



~6 a? 



iJ.O 



120 J? 



ijfi 



0{h% 



(33) 



Assuming sufficient smoothness of v and /, and using OTT l and (|32] |. we can express 
the third and fifth derivatives of v in the form 



53y 

dz^ 
dz^ 



z=0 



z=0 



z=0 



-Aik 



2k^ I- 

dz} dz?dx^ 



z=0 



dz 



2ik 



dz?-dy^ 

dz} dzdx^ dzdy^ dz? 



z=0 



(34) 
(35) 



To preserve the sixth order compact approximation in ( l33l l we need to use the fourth 
order approximation for ^ in ( |34] |. The second order approximation is sufficient for 
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the derivatives of v in ( |35] ). First, let's consider the fourth finite-difference compact 
approximation for It can be presented in the form 



ijfi 



12 



iJ.Q 



12 
12 



4t 



(36) 



2ik 



df d^f 



dz dz?- 



-0{h^ 



2ik 



i.jfl 
dz dz?- 



iJ.O 



0(h% 



Now we can use (l32b . (l36] l and the second order approximation of dSSl) to express 
( l33T l in the form 



ikh^ 



15 

6 ^ 45 / <9z 



^^'^'^ 5>„-,o + ^ (5252v,,,o + 5^5>,-,,,o) =^;-,;,o, (37) 



180 
"90" a? 



ij,0 



120 



^V" d^f d'f 

dz^ dzdx^ dzdy^ 



i.jfi 



This equation gives an implicit compact sixth order approximation for ( |32] ). Next, the 
explicit implementation of the boundary conditions ( |32] ) similar to ( |29] | is developed 
by using the equation 



= /iVi 



dzdx^ dzdy^ 



(38) 



6 \5z^(9x2 d^zdy^ 



i.j.O 



ikjJ-ir \ d'*v d*v 



lixh 



4 r 



2Ai2 



dx^dy^dz dzdx^ dzdy^ 



dz}dx^ d'^zdy^ 
+ 0(h^). 



iJfi 



i.jfi 



Here, /Xj and H2 are parameters. We also use that 



= Oand 



dzdv^ 



- 0. We leave the first of these derivatives in the right hand side but drop the 



latter two derivatives. Now, by using the second order finite-difference approximation 
in ( l38T l and adding it to the (l37l l. we derive the formula 



M3 



OzViJfi - — ( 1 



2l,2 



Ik^h 

~~15 



ikh"^ 



-M3^[5^5^,,o + 5|52v,,,o](39) 



dzdx^ dzdy^ 



iJfi 
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Next we choose parameters /Xi,/i2 and jj.^ to match the coefficients in dZSl l. This 
choice is determined by the conditions 



Ml 



ikh 
90 



1 

6fi3 



kfhf 
30 



1 



2ikh 
3 



1 



10^1^3 



(40) 



Using these parameters and replacing Vi j j with Ui j je -', we write 



„likh 



M3 



/iM Ml + ^ 



30 
2ikh 
3 



(41) 



+ 



+ e'*^V3 



e Us 



30^3 



15 



IhFij.o + liij 



^7 , d\f 



2ikh 4ih^P 
3 

4-ikh^ 



45 
1 



t^.,;.i 



dzdx^ dzdy^ 



60 



i,;,0 



This equation provides an explicit compact sixth-order approximation for the bound- 
ary condition ( l30l l. At the upper boundary z = a, the Sommerfeld-like boundary 
condition can be written ^ — iku — and the auxiliary substitution becomes 

V = e^'^^^u. The rest of the derivation is very similar to the derivation for the case of 
the lower boundary. Implementations of the sixth order compact approximation of 
the boundary conditions in the other directions are also similar to the calculations 
already presented. 



3.3 Fast Fourier Preconditioner 

The main idea of this paper is to utilize the existing lower-order approximation direct 
solvers as preconditioners in the implementation of a higher resolution scheme. One 
of the most popular methods for the approximate solution of the Helmholtz equation 
is the second-order central difference scheme. We consider the cases with Dirich- 
let or Neumann boundary conditions on the sides of the computational domain Q, 

i.e. u{0,y,z) = u{a,y,z) = u{x,0,z) = u{x,a,z) = or |^ 



x=0 



du 



du 

7J- 



du 



= 0. At the top z = a and the bottom z = of the computational domain, we as- 
sign one of the three boundary conditions under consideration: Dirichlet, Neumann or 
Sommerfeld-like boundary conditions In the case of Neumann and Sommerfeld- 
like boundary conditions, we use the second order central finite-difference approx- 
imation of the first derivative on all boundaries. To express the preconditioner in 
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general form we introduce the N x N matrices 

"0 l+j3 



A 



1 1 







1 
l+j3 



(42) 



and D" = diag (a,0,...,0,a), where a , j3 and are parameters depending on bound- 
ai-y conditions. Next let ^ = - 2/ and A^.p^ = A^ - (2 - k^h^jl + D^. The 
preconditioning matrix can now be written in the form 



A„=A 







/ «) / + / (8) A ^ (g) / + / (8) / (Ki A „_ , p. , 



(43) 



where j3 = in the case of Dirichlet boundary conditions at the sides of the compu- 
tational domain and j3 = 1 in the case of Neumann boundary conditions on the same 
boundaries. The number of grid points in the z-direction A^, is chosen in such a way 
that the grid step size is the same in all three directions. This condition is imposed 
so that the lower order preconditioner and the sixth-order approximation compact 
scheme developed in the previous section use the same grid points. Depending on 
the boundary conditions at the top and the bottom of the computational domain, the 
parameters a, and j3, are defined in the following way: in the case of Dirichlet bound- 
ary conditions a, = and jS, = 0; in the case of Neumann boundary condition a, = 
and jSj = 1 ; and in the case of the Sommerfeld-like boundary conditions — 2ikh 
and Pz = 1- At each step of the iterative process, the solution to the second equa- 
tion in (fTTI) with the matrix (|43] | can be obtained by using FFT-type algorithms in 
0{N^NAog2N) operations. 



3.4 Convergence analysis of 3D algorithm 

In this section, the convergence of the proposed algorithms is considered in the case 
of the 3D Helmholtz equation with the Dirichlet boundary conditions ([T]), ^ on 
the rectangular computational domain £2 = {0 < x,y,z < 1} with uniform grid size 
h — ^yq-j- . To insure the uniqueness of the solution to the original boundary value prob- 
lem, we assume that the coefficient is in the resolvent set of the Laplace operator 
defined on an appropriate function space. We also impose a common restriction on 
the number of grid points per wave length, i.e. kh < The convergence analysis of 
the proposed algorithms follows the same lines as in Theorem[3] The preconditioning 
matrix is given by (|43]) with ^ — a, — jS^ — and A^, ~ N. To express the right hand 
side of the sixth order approximation scheme given in ( |25] ), we will use the notation 
introduced in ( |43] ). In addition, we define the N xN matrix A^ ^ = A^ — 47. Then the 
resulting A^^ x A^^ matrix A can be written as 

A ^ A^j^® I ® I + 1 (SA'^j^i^ I + 1 (SAIj^+ (44) 
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1 

6 
1 

30' 



27,2 



h-'k 



(A^ (g) ^ (g) / + / «) A,^ ^ (g) - 4/ (g) / (g) /) - 



12 



360 



/ig)/(g/. 



Both matrices A and Ap have the same set of orthonormal eigenvectors V, j / (see 
e.g.il): 



(45) 



^i")'/'^ ~ i'^h)^^'^ sin (imnh) sin [jnnh) sin (Isnh), 
1 < ijil,ni,n,s < N. 

The matrix of eigenvectors is unitary so the Z2-norm of V and V^^ is one. To present 
the eigenvalues of the matrices, we introduce the notation Sr = sin^ (^),r = 
Then the eigenvalues of the matrices Ap and A in the 3D case are given by 



-4 [.2 



h^k^, and 



-hV- 



(46) 



Ifk^ 
30 



r 2 2 



2 2 



2 21 



32 
15' 



2 2 2 



h^k^ h^k^ 



12 



360 



= 1, ...,A^,respectively. 



Unfortunately, the eigenvalues of these matrices do not satisfy the hypotheses of The- 
orem [T] but we still can apply Corollary [T] to prove the convergence of the GM- 
RES method and the SKS algorithm by showing that Xijj/Xf'j / — 1 +di jj, where 
\dijj \ < 1, for all = 1, ...,A^ in some simple cases. 

Let's assume that > 1, i.e. /i < i andk^ < 25 . Theni^ = sin^{^) > sin^(^) > 



2h , for r — i,j,l. It is also easy to see that ^sfsj > j^sfsjsf and ^{sfsf + sjsj) > 

/■A I 4 /4/4 

2f- > for the given values of k. Now the eigenvalues of the preconditioned 
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matrix AA^ ^ can be written in the form 



3'^ 



30 



■ 2 2 



„2 2 



2 21 



32 
15^ 



2 „2 2 



6;,6 



12 360 



(47) 



This gives us a presentation of the eigenvalues of the preconditioned matrix in the 
form Xijj/Xj'j I = 1 +dij^i, where \di j i\ < 3/4, for all i,j,l = l,...,N. Now we can 
apply Corollary [Hand obtain an estimate for the convergence rate of convergence of 
the GMRES and SKS algorithms : 



< (3/4)"||ro||2, 



(48) 



where r„^F-AU^"'. 

Similarly to the ID case, it is possible to show that this estimate holds for ^ > 6 
and a sequence of grids with the grid sizes hi > h2 > . . . if 



mm 

i,jJ=l.....,N 



sj + sj + sj 



/l2 



>5o>0 



(49) 
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for some 5o — const and sufficiently small hi. As in the case of the ID problem if 
follows from the fact that the resolvent set of a bounded linear operator is open. The 
proof is somewhat tedious but elementary, so we skip it. In the next table we present 
lower and upper bounds m and M for A, j./ /Xfj ^ — 1 for different values of k and for 
different grid step sizes. Also, the last two rows of the table present possible choices 
for hi and 5q for different k. 

Table 1. 3-D Dirichlet problem eigenvalue bounds(different k^) 





k^lQ 


A: = 20 


k = 3Q 


;t = 40 


A: = 50 


h=\l6A 


-.015 10.49 


-1.22 1 .49 


-1.02 1 .88 


-34.7 1 180 


-10.2 1 55.0 


/z= 1/128 


-.0038 1 0.49 


-.158 1 .49 


-.140 1 .64 


-2.81 1 0.90 


-2.21 1 1.74 


h = 1 /256 


-.00094 1 0.49 


-.035 1 .49 


-.031 1 .49 


-0.22 1 0.69 


-11.7 1 .56 


/i= 1/512 


-.00023 1 0.49 


-.009 1 .49 


-.008 1 .49 


-.048 1 0.49 


-.298 1 .49 


do 


8.00 


4.00 


1.80 


1.12 


2.9 


hi 


1/64 


1/128 


1/256 


1/256 


1/512 



This result is weaker then in ID case, i.e. it says that convergence is independent 
of the grid step, but it does not improve with the decrease of h. In this case, Ap is said 
to be "an optimal preconditioner" (see e.g. lil5J . p. 196). This situation is typical for 
multigrid-type algorithms but it is not what we would expect based on \D example. 
In the following section we present the results of numerical experiments in the case of 
several test problems which confirm that actual convergence of the presented methods 
accelerate with the decrease of the grid size, i.e. the convergence in the numerical 
experiments significantly exceeds the estimate ( |48] |. 

4 Numerical Results 

In this section we present the results of numerical experiments which demonstrate 
the quality of the proposed numerical methods. These algorithms were implemented 
in Matlab 7. 1 1 .0 on an iMac with an Intel Core 17, 2.93 GHz processor. We also use 
the standard programing implementation of the GMRES method in Matlab (gmres 
function). 

4.1 ID test problem 

In the first series of numerical experiments, we consider the convergence of the 
compact sixth order approximation algorithm on a sequence of grids. In these tests 
we focus on the numerical solution of the ID Dirichlet problem for the Helmholtz 
equation with zero boundary conditions on the interval < x < 1 and k = 20. The 
computational grid in our experiments is defined by x, = //;,/ = 0, ...,A^ + 1 where 
h = ^yJpY ™d is the number of grid points. We consider the function m(x, ) = x, (1 — 
Xi)cos{knxi),i = 1,...,N as the exact solution of the original boundary value prob- 
lem with the right hand side /(x,) = —{2 + k^{7i^ — x, (l — x,)) cos(A:;rx,) +2A;;r(2x, — 
l)sm{k7ZXi),i = l,...,N. The first column of Table 2 displays the step size h of the 
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grids in our experiments. Columns 2-4 of the table present the number of iterations 
until convergence required by the preconditioned GMRES method, the SKS method 
(Algorithm[T]i and the Chebyshev acceleration (CA) algorithm [ 19 |. We use the stop- 
page criterion ||F — At/("'||2/||f'||2 < 10^'", where t/'"' is the iterative solution on 
nth iteration. Column 5 reports the approximation error of the numerical solution 
Errs = maxi<j<A? jt/j"' — u{xj)\. In Column 6, the second order approximation er- 
ror Err2 of the numerical solution obtained by using the preconditioner as a solver 
is presented for comparison. The number of grid points is chosen to satisfy the re- 
quirement h < 2n/lQk. The values of the minimal and maximum eigenvalues of the 
preconditioned matrix used by the Chebyshev acceleration algorithm are calculated 
exactly by ( l22b . 



Table 2. 1-D Dirichlet problem 



Grid Size (h) 


GMRES 


SKS 


CA 


Errs 


Err2 


1/32 


6 


8 


6 


3.88*10-^^ 


1.16*10-*'^ 


1/64 


4 


4 


5 


5.08*10-"^ 


2.4*10-"^ 


1/128 


3 


3 


4 


7.61* 10-' 


5.7*10-"^' 


1/258 


3 


3 


3 


1.14*10-** 


1.4*10-"^ 



These experiments demonstrate the sixth order convergence of the numerical so- 
lution to the exact solution of the original boundary value problem. Theorem |3] sug- 
gests that the residual on nth iteration in the presented iterative algorithms decreases 
with the increase in the number of grid points. The results of the numerical experi- 
ments confirm this conclusion. We have already shown that the order of the precon- 
ditioned matrix in this case is two. The numerical value of the constant y/ calculated 
by using (l23l l is 1.99. The ID numerical experiments confirm the effectiveness of 
the iterative approaches under consideration and allows us to expect that the same 
properties hold for 3D problems. 



4.2 3D test problems 

In the next series of experiments we consider different boundary value problems 
for the 3D Helmholtz equation ([T]!,©. The computational domain is given hy £2 — 
{{x,y,z)\Q < x,y,z < 1 } with a uniform grid in all three directions Xi — ih,i — I, ...,Nx, 
yj = jh,j = l,...,Ny andz^. = kh,k= l,...,N^. The numbers of grid points A^^jA'j and 
in x,y and z-directions are chosen to provide uniform grid size. 

4.2.1 Dirichlet problem 

First, we consider the Helmholtz equation with boundary conditions m = and k = 20. 
This test is similar to the ID case with A^^ = Ny = N, = N and h = 1/(A^+ 1). 
As the solution of the boundary value problem we choose the function u{x,y,z) = 
01 {x) 02 (y) 03 (z) , where (j)i{x)=x^{l-x)^,(j)2iy)=y{l-y)cos {kny) , (z) = sin(^7rz) , 
and < x,y,z < 1. 
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In the first series of tests we consider the convergence of the compact high order 
approximation finite-difference scheme dZSl l on a sequence of grids as well as con- 
vergence properties of the GMRES, SKS and CA iterative algorithms. The iterative 
process is stopped when the initial residual is reduced by a factor of 10^ The values 
of the minimal and maximumal eigenvalues of the preconditioned matrix used by the 
Chebyshev acceleration algorithm are calculated exactly using ( |47] |. The description 
of the columns of Table 3 is similar to the description of Table 2. 



Table 3. 3-D Dirichlet problem 



Grid Size (h) 


GMRES 


SKS 


CA 


Erre 


Err2 


1/64 


5 


10 


23 


3.03 * 10-*"^ 


4.47 * 10-** 


1/128 


3 


6 


15 


8.13*10"''^ 


6.35 * 10-*'** 


1/256 


3 


5 


14 


2.06*10-''^ 


9.68*10-'** 



The test results confirm the sixth order approximation of the compact finite-difference 
scheme. Though, in our theoretical analysis ( |47] | we could not prove the desired kth 
order preconditioning property in the case of the 3D Dirichlet problem, the SKS 
algorithm still exhibits acceleration of the convergence with the decrease of grid size. 
We also mention that the total computer time required for convergence for the SKS 
and GMRES methods is approximately the same and the Chebyshev acceleration 
algorithm requires about four times more computer time until convergence than the 
other two iterative algorithms do. The SKS algorithm has also greater potential for 
efficient implementation on parallel computers than the GMRES method does. 

One potential application of the proposed methods is to electromagnetic scat- 
tering problems. In these problems, the typical value of the coefficient k is in the 
range between 10 and 50. For example, the propagation of an electromagnetic wave 
with 1 GHz frequency in a vacuum is described by the Helmholtz equation with 
sa 20.9 m-'. In the next table we present the convergence of the algorithms pre- 
sented here on the sequence of grids with different values k. In every cell of the table 
we present the number of iteration for three methods under consideration ; GMRES, 
SKS, and CA. The stoppage criterion is the same as in the previous series of exper- 
iments. If any of these three methods fails to converge in 100 iterations, we indicate 
this by using symbol "> 100" and in the case of divergence we use "div" symbol. 

Table 4. 3-D Dirichlet problem (different k^) 





k=\Q 


k = 2Q 


k = 3Q 


k = AQ 


;t = 50 


h = 1 /64 


4| 7 1 14 


5| 10 1 23 


6| 13 1 25 


12| div 1 > 100 


19| div 1 > 100 


/i= 1/128 


3| 5 1 14 


3| 6 1 15 


4| 8 1 17 


4| 10 1 69 


5| 11 1 > 100 


h = 1 /256 


3| 4 1 12 


3| 5 1 14 


3| 6 1 13 


3| 6 1 22 


4| 7 1 63 




2.02 


2.00 


1.98 


2.00 


2.02 



From these numerical experiments we can see that the preconditioned GMRES and 
SKS methods demonstrate excellent convergence properties. We also observe that the 
convergence of the proposed algorithms improves with the increase in the number of 
the grid points. In the analysis of the algorithm's convergence (l48T l. it was shown 
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that the number of iterations does not increase with the increase of the grid points. 
But the last row in the table indicates that the parameter \(f ( l23T l for all values of k 
is approximately 2, which indicates that similarly to the ID case the convergence of 
the 3D algorithm not only does not depend on the grid size but accelerates with the 
increase of the number of grid points. This also indicates that even in the case when 
the SKS method does not converge, one can expect convergence on a finer grid. In all 
numerical tests, the parameter y/ is calculated by using the first two iterations of the 
SKS method on the grids with the grid sizes h = and h = • 

From the results presented, it is clear that slow convergence and difficulties arising 
in the choice of the parameters m and M in the Chebyshev acceleration algorithm (fTSl ) 
make this method a poor choice for the implementation of the developed high-order 
approximation scheme. These numerical tests illustrate the fact that the Chebyshev 
polynomial is not optimal on the discrete spectrum and this has a dramatic effect on 
the iterative methods based on this polynomial. So, in the next numerical experiments 
we focus on the convergence properties of the GMRES and SKS methods in our 
numerical framework. 

4.2.2 Dirichlet-Neumann boundary conditions 

In the next series of experiments we replace the Dirichlet boundary conditions u — Q 
in (|2| with u — Uj, at the top of the computational domain (z = 1) and with the 
Neumann boundary conditions ^7 = at z = 0. On the other boundaries of the 
rectangular computational domain £2 = {0 < x,y,z < 1} we still use the Dirichlet 
boundary conditions m = 0. As a test function that satisfies the boundary conditions 
we consider u{x,y,z) = ^i{x)^2(y)^3iz), where =x^(l —x)'^,^{y) — y{l — 

y)cos{kny),^{z) — cos{k7rz), and < x,y,z < 1. Also, Uj, — ^i{x)^{y)^^{l),0 < 
x,y < 1 . The number of grid points is chosen to provide the same grid step in all three 
directions, i. e. Nx = Ny = — The stoppage criterion is the same as in the previ- 
ous numerical experiments. Table 5 provides the results of the numerical tests on a 
sequence of grids with the coefficient k ~ 20. As already mentioned we focus only 
on two best algorithms: GMRES and SKS. The rest of Table 5 is similar to Table 3. 



Table 5. 3-D Dirichlet-Neumann problem 



Grid Size (h) 


GMRES 


SKS 


Erre 


Err2 


1/64 


8 


11 


4.65 * 10-'"' 


3.67 * 10-'''* 


1/128 


5 


8 


6.61 * 10-*'** 


8.74*10-''^ 


1/256 


4 


7 


1.01 H^IO-*'^ 


2.17*10-''^ 



The results of this series of numerical tests indicate the sixth order convergence of the 
resulting compact finite-difference scheme including explicit approximation of the 
Neumann boundary conditions. As expected, the number of iterations decreases with 
the increase of the number of grid points in both the GMRES and SKS approaches. In 
the next table we present the number of iterations required for convergence for both 
algorithms on a sequence of grids and with different k coefficients in the Helmholtz 
equation. The content of Table 6 is similar to the content of Table 4, except for the 
number of iterative approaches under consideration. We present only the number of 
iterations for the GMRES and SKS methods. 
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Table 6. 3-D Dirichlet-Neumann problem (different k^) 





k=\Q 


yt = 20 


k = 2>Q 


;t = 40 


k = 5Q 


h=\/6A 


6 1 10 


8 1 11 


18 1 div 


36 1 div 


>600 1 div 


h= 1/128 


4 1 9 


5 1 8 


8 1 9 


14 1 div 


28 1 div 


h= 1/256 


3 1 8 


4 1 7 


5 |div 


9 1 div 


7 1 div 


1/512 


2 1 6 


3 1 6 


3 1 6 


5 1 6 


4 1 6 


W 


2.57 


2.18 


2.39 


1.65 


1.94 



It follows from this table that the GMRES approach is much more efficient on coarse 
grids but on finer grids, the processor time for the SKS method is the same or even 
smaller than for the GMRES algorithm. For example, on the grid with h = 1/512 
and k = 5Q the corresponding processor times for the GMRES and SKS methods are 
1467 sec and 1052 sec. Also, as we mentioned before, the SKS method has much 
greater potential than the GMRES algorithm for efficient implementation on parallel 
computers. We also expect that even if the SKS method does not converge on coarser 
grids, by reducing the grid size one can achieve convergence of this algorithm. We 
can see this in the example of the convergence of the methods under consideration 
in the case of A: > 30. The values of the parameter vary significantly for different 
values of k but > still indicates that the convergence of the iterative algorithms 
improves with the increase in the number of the grid points. 

4.2.3 Dirichlet-Sommerfeld-like boundary conditions 

In the last series of numerical experiments, we consider the boundary value problem 
for the Helmholtz equation ([T]) with a combination of Sommerfeld-like(radiation) 
boundary conditions (|2]l at z = and z = 1 , and Dirichlet boundary conditions at all 
other boundaries of the rectangular computational domain £1 ^ {Q < x^y^z< 1}. The 
source function / in ([T]) selected that the true solution is u{x,y,z) — 0i {x)^2{y)^i{z), 
where ^ x^{\ - xf ,<^2{y) ^ y{l ~ y)cos{k7Zy),(j)3{z) = e'<^('+i) +e"'^'(-'"') -2 
and < x,y,z < I. Note that the analytic solution satisfies the radiation boundary 
condition (|2]). The numbers of grid points in XjjjZ— directions are chosen such that 
A^T = A^v = N^ — l. First we consider the convergence of our sixth order approxi- 
mation scheme on the sequence of grids with k = 20. The main goal of this series of 
numerical experiments is to investigate the convergence properties of the new explicit 
compact sixth order approximation of the Sommerfeld-like boundary conditions (HTI) 
proposed in this paper As in previous test runs, the iterative process is stopped when 
the norm of the initial residual is reduced by a factor of 10^^". The results are 
presented in Table 7. The description of data presented in this table is the same as in 
the case of Table 5. 



Table 7. 3-D Dirichlet-Sommerfeld problem 



Grid Size (h) 


GMRES 


SKS 


Errs 


Err2 


1/64 


6 


9 


1.25*10-*'^ 


1.35*10-''^ 


1/128 


4 


7 


1.84*10-*'** 


3.23*10-""* 


1/256 


3 


5 


2.82*10-'" 


7.99*10-"^ 
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In this series of numerical experiments, both algorithms exhibit the same convergence 
properties as in the previous series of tests, i.e. the number of iterations decreases as 
the number of grid points increases. From the table, we can also see the sixth order 
convergence of the approximate solution to the exact solution of the boundary value 
problem on a sequence of grids. These results confirm the sixth order approximation 
of the compact explicit scheme proposed for the numerical implementation of the 
Sommerfeld-like boundary conditions. 

In the last table, we present the number of iterations required for the convergence 
of the GMRES and SKS methods for different coefficients k in the Helmholtz equa- 
tion with the same boundary conditions used in the previous series of numerical tests. 
Data presented in the next table are sinoilar to the data in Table 6. 

Table 8. 3-D Dirichlet-Sommerfeld problem (different A^) 





k = 10 


k = 20 


k = 30 


k = 40 


k = 50 


k = 35.7 + 0.43/ 


ft = 1/64 


5|7 


6 1 9 


11 1 43 


16 1 div 


119 1 div 


14 1 54 


ft = 1/128 


4|5 


4 1 7 


7 1 8 


7 1 div 


12 1 11 


7 1 9 


ft =1/256 


3 1 4 


3 |5 


5 1 6 


5 1 6 


6|7 


5 1 6 


ft= 1/512 


2 1 4 


3|5 


4| 5 


3 1 5 


4 1 6 


4| 5 


¥ 


15 


13 


13 


15 


11 


13 



In the last column of Table 8, we present the results of calculations when k is a 
complex valued coefficient. The value of this coefficient approximately corresponds 
to the case of propagation of electromagnetic waves with the 1 GHz frequency in dry 
soil. These results confirm the effective implementation of the compact sixth order 
approximation scheme by using the proposed Krylov subspace-type algorithms in the 
framework with the FFT based low order preconditioners. The last row of the table 
suggests that the acceleration of the convergence with the increase of the number of 
grid point is much stronger than the acceleration in the previous series of experiments. 
It seems strongly dependent on the boundary conditions used in the numerical tests. 
In all our experiments, there is a clear connection between the parameter \j/ and the 
improvement of convergence with the decrease of the step side. But the usefulness 
of this parameter in the analysis of the quality of a preconditioner requires further 
analysis since in majority of our numerical experiments the preconditioner does not 
satisfy the condition in the Definition 1 . 

The series of test problems considered suggests that in the majority of situations 
the preconditioned GMRES method is the most efficient choice for an effective im- 
plementation of the compact sixth order approximation scheme on the coarse grids 
but in the case of finer grids the SKS method in combination with lower order approx- 
imation preconditioner presents an efficient alternative to the GMRES method. This 
alternative could become even more valuable when the GMRES method experiences 
stagnation. 
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5 Conclusions 

New 3D compact sixth order explicit finite-difference schemes for the approximation 
of Neumann and Sommerfeld-Uke boundary conditions on rectangular computational 
domains with uniform grid size were developed and implemented. Together with the 
compact sixth order approximation scheme for the Helmholtz equation proposed in 
En . these algorithms represent highly accurate methods for the solution of boundary 
value problems for Helmholtz equations. 

A new rapid iterative method based on preconditioned Krylov subspace method- 
ology was developed for the implementation of the proposed compact finite-difference 
schemes. The strategy is based on a combination of higher order approximation 
schemes and a lower order approximation preconditioner The analysis of some typ- 
ical test problems reveals the attractive properties of the developed methods such as 
the decrease of the number of iteration until convergence with the increase of the 
number of the grid points or the size of the resulting matrix. This approach is espe- 
cially attractive in situations in which the lower approximation solver already exists 
and the original boundary value problem calls for more accurate approximation. 

The typical time to produce the sixth order accuracy solution of the 3D Helmholtz 
equation with a combination of the Dirichlet and Sommerfeld boundary conditions 
on a 512^ grid was just 30 minutes on a iMac using only one 2.93 GHz Intel Core 
i7 processor. The method was tested for realistic parameter ranges typical for elec- 
tromagnetic scattering problems. We must notice that the Sommerfeld-like boundary 
is just a first order approximation of the Sommerfeld conditions on the unbounded 
domain. So, direct application of the higher order approximation for the Sommerfeld- 
like boundary condition to the solution of a scattering problem is not always justified. 
However, a straightforward extension of the approximation approach presented in 
this paper could be applied for approximation of the absorbing boundary conditions 
(see e.g. |23 |) or can be used in the implementation of the perfectly matched layer 
(PML) boundary conditions (see e.g. I.13J ). 
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